The semiflexible fully-packed loop model and interacting rhombus tilings 
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Motivated by a recent adsorption experiment [M.O. Blunt et al, Science 322, 1077 (2008)], 
we study tilings of the plane with three different types of rhombi. An interaction disfavors pairs of 
adjacent rhombi of the same type. This is shown to be a special case of a model of fully-packed loops 
with interactions between monomers at distance two along a loop. We solve the latter model using 
Coulomb gas techniques and show that its critical exponents vary continuously with the interaction 
strenght. At low temperature it undergoes a Kosterlitz-Thouless transition to an ordered phase, 
which is predicted from numerics to occur at a temperature T ~ 1WK in the experiments. 



Byzantine artists were probably the first to appreci- 
ate the beauty of tiling the plane with various types of 
polygons. Nowadays, random tilings arc a major sub- 
ject in mathematics and theoretical physics, with many 
applications to real condensed matter systems. 

In a very recent experiment, Blunt et al. [l[ study the 
adsorption of certain rod-like organic molecules (para- 
terphcnyl-3,5,3',5'-tetracarboxylic acid) on a graphite 
substrate by scanning tunneling microscopy (STM). An 
idealization of the resulting pattern is shown in Fig. [T]A 
The rod-like molecules (shown as thick white lines) ar- 
range as a dimer covering of the hexagonal lattice, and 
interact with their neighbors via hydrogen bonding of 
carboxylic acid groups (thick black lines). The bisectors 
of the latter define a rhombus tiling with three different 
types of rhombi. By image processing the STM pictures, 
large tiling configurations could be produced [lj, with 
only ~ 10~ 3 defects per molecule at room temperature. 

Dimer coverings of various planar lattices have been 
well studied theoretically and the corresponding par- 
tition functions and various correlation functions have 




FIG. 1: (color online) Dimer covering of the hexagonal lat- 
tice (thick white lines), the complementary fully-packed loops 
(thick black lines), and corresponding rhombus tiling of the 
plane (A). Height mapping for the rhombus tilings (B) and 
for the fully-packed loop model (C). 



been obtained exactly 0, 0] • Long-distance behavior is 
conveniently analysed in terms of an equivalent height 
model. For the hexagonal-lattice case, this is obtained 
Q by defining a height difference ±1 to each displace- 
ment along the junction between two rhombi, using the 
rule in Fig.[Tj3. In the continuum limit, one expects [3|, Q 
the height to behave as a Gaussian free field with some 
coupling g which can be computed exactly [3, Q . 

The experimentally measured value of g [l| is however 
1.66(8) times larger than the theoretical prediction. This 
discrepancy is attributed to interactions between pairs of 
neighboring dimers, with an energy penalty AE > for 
a parallel arrangement. This motivates the study of in- 
teracting rhombus tilings (IRT), where each tile junction 
carries a weight w = c~ AE / T (respectively 1) if it sepa- 
rates rhombi of identical (resp. different) types. 

The purpose of this Letter is to study a more general 
model which contains the IRT as a special case. We de- 
fine a link to be any edge of the hexagonal lattice not 
covered by a dimer. In Fig. [T]A_ links are shown as thick 
black lines. With appropriate boundary conditions, the 
links form fully-packed loops, i.e. loops which jointly visit 
each of the lattice vertices once. Assign a weight n to each 
loop. Since each link Cq is a bisector of a tile junction, 
the weight w must be attributed to Cq if and only if the 
two links touching either end of Cq are parallel. This 
defines the partition function of the semiflexible fully- 
packed loop (SFPL) model: 
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When n = 1 we recover the IRT. 

Below we solve the SFPL model in the critical regime 
— 2 < n < 2 using Coulomb gas methods 0, Q. The 
case w = 1 was previously shown [j| to be described 
by two free fields. We shall see that the interaction w 
changes the coupling g\ of one of these fields, leaving 
the other gi unchanged. The ratio 7 = g\jgi is a non- 
universal decreasing function of w. We compute all criti- 
cal exponents as functions of n and 7. When w becomes 
smaller than some critical w c , the SFPL model under- 
goes a Kosterlitz-Thouless (KT) transition to an ordered, 



non-critical state. Below we compute analytically the 
corresponding value of j c , and perform numerical simu- 
lations to compute the curve j(w). This allows to predict 
the value of the temperature at which the KT transition 
should occur in the experiments. We conclude the Letter 
by a critical discussion of the experimental realization [l| , 
and a comparison with other related loop models. 

2D height mapping. We briefly review the construction 
of Refs. p, 0] and show how elements must be modified 
to account for the interaction w in the SFPL model. 

Orient each loop independently and assign a weight 
e -i7reo £ clockwise and e I7re ° to counterclockwise loops. 
Parametrizing n = 2cos(-7reo) we recover the correct loop 
weight after summing over orientations. Note that this is 
tantamount to a weight e ±17re °/ 6 to each left (resp. right) 
turn. Assign a label vq to edges covered by a dimcr, 
and v± to edges covered by a link going out of (resp. 
into) a vertex in the even sublatticc. Each vertex is 
then adjacent to three edges, all carrying different labels 
(vq,v+,v-). Define the corresponding two-dimensional 
vectors v = (—2,0) and v± = (1,±VS) (see Fig. [Tp). 
Attribute 2D heights h = (h 1 , h 2 ) to the dual triangular 
lattice, by increasing h by v, upon traversing an edge 
with label i. The traversal must be such that an even 
(resp. odd) vertex is seen on one's left (resp. right). The 
first component h 1 is precisely the ID height defined by 
Fig. [TJ3 for the rhombus tiling. Being complementary to 
loops with no orientation, the rhombi cannot "see" h 2 . 

Coulomb gas approach. The partition function is writ- 
ten as a functional integral 

Z = J Ph(x) cxp(-S[h(x)]) , (2) 

where, by an abuse of notation, h(x) denotes the contin- 
uum limit of the height defined above, and the Euclid- 
ian action consists of three terms, S = Se + Sb + Se- 
The elastic term Se is constrained by rotational invari- 
ance to take the form [with summation over repeated in- 
dices] S E = 1/2 fd 2 yig a pdh a ■ dh , where d = (d 1 ,d 2 ) 
is the usual gradient. The D = 2-dimensional sym- 
metric tensor g a p is further constrained by symmetries. 
First, since loop orientations arc eventually summed over, 
the action must be invariant under v+ — > V—, viz., 
{h},h 2 ) — > (h 1 ,— h 2 ), implying g± 2 = 0. We denote 
henceforth g\ = g\\ and 172 = 322- Second, a cyclic per- 
mutation of (vq , U+, u_) maintains the chirality of the 
loop turns at each vertex, and is thus a symmetry for 
w = 1. This implies g\ = g 2 in this case 

However, the cyclic permutation is not a symmetry of 
the SFPL model for w ^ 1. To see this, we inspect all 
possible configurations of the pair of vertices surround- 
ing a fixed edge £q. By rotation symmetry we can take 
£0 horizontal with an even vertex on its left end. By re- 
flection symmetry it suffices to inspect six out of twelve 
configurations. This gives two groups of three configu- 
rations related by the cyclic permutation. As seen from 
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FIG. 2: The upper group of three configurations for a vertex 
pair has weights \ = e l7Te °^ 3 which are invariant under the 
cyclic permutation of labels (vo,v+,V-), whereas the weights 
of the lower group differ by factors of w. 




FIG. 3: Each column shows a pair of configurations with iden- 
tical link positions, but different loop orientations. The mid- 
dle column comes with a weight w. 

Fig. [2] the members of the second group differ by fac- 
tors of w, proving the statement. We now argue that 
changing the weight w will modify g\ but leave g 2 un- 
changed. To this end, we consider pairs of local config- 
urations having identical link positions, but which differ 
by the loop orientation. In Fig. [3] we show three such 
pairs (all others can be found by reflection and rotation) 
and evaluate the height gradient along the middle edge. 
In all cases, the w interaction must not distinguish be- 
tween two members of a pair. We conclude that w must 
not couple to v + — v = (0, 2\/3), which is proportional 
to h 2 . Conversely, configurations in the middle column 
of Fig. [3] have weight w and are the only ones to have a 
height gradient along Vo cx h . This proves the claim. 

The action also contains a boundary term Sb = 
i/(47r) J d 2 x (eo ■ h)72.(x) where TZ is the scalar curvature. 
The background electric charge eo is easily computed 
on the cylinder, where it ensures the correct weighting 
of non-contractiblc loops. We have eo ■ vo = and 
e • V± = ±7re , implying e = (0,7re /V3). 

Finally, the Liouville term is the continuum limit 
of the local vertex weights. The height is compactificd 
[j| with respect to a triangular lattice M of side 2\/Z 
spanned by v± — v . We can therefore expand as a 
Fourier series over the vertex operators e le h , where the 
electric charges e belong to the lattice £ reciprocal to M. . 
It suffices to keep the most relevant term which has the 
same periodicity as the vertex weights. The inclusion of 
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w does not change this periodicity, so we have Sx ~ e I6s ' h 
where the screening charge reads e s = (0, 2n/y/3). 

Critical exponents. The critical exponent of an opera- 
tor O e m with electric charge e £ £ and magnetic charge 
me M reads [1, 0] : 

x e , m = [e a (e a - 2e Qa )/g a + g a m 2 a ] /(An) . (3) 

The corresponding two-point functions decay with dis- 
tance r as r~ 2x °- m . To keep the model critical we impose 
Q the exact marginality of S*l, whence x es ,o = 2, or 

0j = (l-e o )7r/6. (4) 

The other coupling g\ depends non-universally on the 
microscopic weight w. Henceforth, we express everything 
in terms of eo and the ratio 7 = 31/(72- Note that 7(11?) is 
monotonically decreasing, and 7(1) = 1. Computing the 
actual function j(w) would require an exact solution, but 
it is doubtful that the SFPL model is integrable. Below, 
we perform numerical simulations to obtain this curve in 
the n = 1 case relevant to experiments. 

The central charge is determined from the background 
electric charge as c = 2 + 12a- en . = 2 — 6eg/(l — eo), 
independent of 7. 

An important class of critical exponents is the so-called 
watermelon exponents Xk- They measure the probability 
of having k oriented loop strands emanating from some 
small neighborhood (of size a few lattice constants) and 
absorbed by some other neighborhood at distance r>l. 
The corresponding height defect (vortex) has magnetic 
charge which is computed by noting that v + — v 
generates a pair of strands, and 2vo — v + generates a 
single strand [f|. Setting 4 = k mod 2 G {0, 1} we find 
explicitly = (— 34, s/3k). This leads to 

x k = x eo ,m k = {k 2 + 3 7 4) (l-eo)/8-eg/(2-2eo) • (5) 

Another type of vortex corresponds to having a vertex 
not visited by any loop. The corresponding magnetic 
charge is m<r = 3vo = (—6, 0), and the exponent 

xt = xo,m T = 37(1 - e )/2 . (6) 

In the IRT model, the second height component is 
"invisible" and only x\ and xt are meaningful. In the 
tiling picture, the corresponding defects are compounds 
of three (resp. four) elementary triangles forming a trape- 
zoid (rcsp. triangle) of base length two. Since eo = 1/3, 
we have gi — tt/9, c = 1, xi = 7/4 and xt = 7- 

Kosterlitz-Thouless transition. Due to the compacti- 
fication, any functional of the heights can be expanded 
over vertex operators with charges in £, a triangular lat- 
tice of side 27r/3. The crucial step in solving the SFPL 
model was to fix the coupling by requiring the exact 
marginality of 5l- By examining the symmetries of local 
vertex weights, one finds Q that vertex operators appear- 
ing in the expansion of Sl have charges in a sublatticc 



£l C £, a triangular lattice of side 2tt/^/3 spanned by 
the second-shortest vectors in £ . For eo > 0, the most 
relevant vertex operator has e = e s = (0,2n/y/3), and 
the solution Eq. ^ was obtained by using this as the 
screening charge, i.e. by setting a; esi o = 2. 

One can add another term Sf to the action that favors 
domains where the height interface is locally flat. Its ver- 
tex operators have charges in a sublattice £f C £l, a tri- 
angular lattice of side 2-k spanned by the second-shortest 
vectors in £l- The interface is in a rough, critical (resp. 
smooth, non-critical) phase when Sf is irrelevant (resp. 
relevant). Increasing 7 beyond a certain critical value 
7 C > 1 induces a KT transition to the smooth phase. The 
most relevant vertex operator in Sp has ep = (±27r, 0). 
We can find j c by setting x eF ,o = 2, yielding 

7c = (l-eo)/3. (7) 

In the non-critical phase, the couplings g\ and gi will 
renormalize to infinity, corresponding to a microscopic 
parameter w — > 0. The situation w = corresponds in 
Fig. UK to the links forming small loops of length six 
around one of the three sublattices (denoted C a with 
a = 1, 2, 3) of the triangular lattice. Which sublattice is 
selected is a matter of spontaneous symmetry breaking. 
Equivalently, in the rhombus picture, the average num- 
ber of rhombi 7V a touching a vertex of C a will saturate 
to (Ni, N2, N3) = (6, 3, 3) or any permutation thereof. 

The exact exponents at the KT transition are found by 
inserting Eq. ([7]) into Eqs. §5§ and ©. Note in particular 
that St (7c) —9/2 is independent of eo. 

Experimental realization. As mentioned in the intro- 
duction, a handsome experimental realization of the IRT 
model appeared recently [l| . The energy scales of the ex- 
periment are such that, once formed, the tilings are static 
(but certain defects can move around dynamically; see 
below) . Statistics on the height fluctuations can however 
be obtained by taking STM pictures of various regions. 
These fluctuations were reported to be critical [H, and 
an experimental value 7 cxp = 1.66 ± 0.08 was observed. 
Since 7 cxp < 7 C with j c = 9/2 in the eo = 1/3 IRT case, 
our analysis confirms that the system is indeed critical. 

To fix the experimental energy scales within the IRT 
model, we perform numerical simulations with Transfer 
Matrix (TM) and Monte Carlo (MC) techniques similar 
to those developped in Ref. Q. Fig. shows c as a 
function of w, as determined from the TM calculations. 
A clear c = 1 plateau appears, corresponding to the crit- 
ical phase for w > w c . We also measured 7 from the 
determination of X\ in TM simulations and winding num- 
ber fluctuations in MC simulations [1, Q . The resulting 
7(u>) curve is displayed in Fig. [4j3 . Extrapolations of w c 
to the thermodynamic limit are made: (A) from the TM 
data by studying the intersections c = 1 and 7 C = 9/2, 
giving w c = 0.635(2) in both cases, and (B) from the 
MC data by styding order parameter fluctuations, giv- 
ing w c — 0.640(5). We also find that the experimental 
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FIG. 4: (color online) Numerical simulations of the IRT 
model. As functions of w. (A) Central charge c, obtained from 
three-point fits of the free energy on cylinders of circumfer- 
ence L. (B) Coupling constant ratio 7, found from two-point 
fits for the critical exponent x\ (TM data), and from winding 
number fluctuations on a 768 x 256 sample (MC data). The 
long-dashed vertical line denotes w c = 0.635. 



value 7c Xp = 1.66(8) corresponds to w oxp = 0.845(15), al- 
lowing to determine the energy scale of nearest-neighbor 
interactions as AE ~ 4.25meV. We therefore predict the 
KT transition to occur at T = T c ~ 110-fr' in the ex- 
perimental compound. The transition could be observed 
by monitoring 7 oxp up to the temperature where it takes 
the value 7 C , as in Ref. [l(J For the precise compound of 
Ref. this will require performing the experiment in 
vacuum, to avoid that the solvent freezes fill ] - 

Among the two possible topological defects mi and 
my in the IRT model, the former is by far the most 
probable, since xt = 4a; 1 (=7). Defects of the my type, 
if observable, would indeed be very closely bound. In the 
dimer language, the mi defect can correspond either to 
zero or two dimers incident to the same vertex. Both 
possibilities were observed in the STM scans, although 
the latter was dismissed as a transient image artifact (see 
in particular Fig. 3E of Ref. [l[). The dynamics of defect 
pairs should make it experimentally possible to gather 
statistics on their relative separation r. Given 7 CX p, the 
above theory predicts the corresponding power law. 

Discussion. We have solved a model of semiflexible 
fully-packed loops, and shown that the bending rigid- 
ity couples to just one of the two coupling constants in 
the equivalent 2D height model. Although wc have here 
given a microscopic argument that only <?i was affected, 
it should be noted that this is also a consequence of the 



field theory. Indeed, since the screening charge e s is in 
the h 2 direction, 52 is in fact bound to renormalize to 
the universal value Eq. Q. The field theory should re- 
main valid for other microscopic interactions that have 
the effect of rendering the height interface stiffer. 

The particular case of interacting random tilings ob- 
tained when n — 1 has a physics similar to that of dimer 
coverings of the square lattice with local aligning inter- 
actions 0. The SFPL model can also be compared to 
the 3D height construction used in Ref. to solve the 
Flory model of protein melting. 

Adding a finite density of mr type defects to the SFPL 
model induces a flow towards the well-known dense phase 
of the 0(n) model [l3j]. Starting from the 2D height map- 
ping, the h 1 component now becomes massive, and the 
dense phase is described by h 2 . Since the corresponding 
coupling 52 is insensitive to w, we deduce that bending 
rigidity is irrelevant in the dense 0(n) model and merely 
renormalizes the effective monomer length. 

Note added. Another experimental realization of the 
IRT model has appeared very recently 
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